The generalized contact process with n absorbing states 
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We investigate the critical properties of a one dimensional stochastic lattice model with n (permu- 
tation symmetric) absorbing states. We analyze the cases with n < 4 by means of the non-hermitian 
density matrix renormalization group. For n — 1 and n = 2 we find that the model is respectively in 
the directed percolation and parity conserving universality class, consistent with previous studies. 
For n = 3 and n = 4, the model is in the active phase in the whole parameter space and the critical 
point is shifted to the limit of one infinite reaction rate. We show that in this limit the dynamics 
of the model can be mapped onto that of a zero temperature n-state Potts model. On the basis of 
our numerical and analytical results we conjecture that the model is in the same universality class 
for all n > 3 with exponents z = «^||/«^± = 2, i^x = 1 and (i = 1. These exponents coincide with 
those of the multispecies (bosonic) branching annihilating random walks. For n = 3 we also show 
that, upon breaking the symmetry to a lower one {Z2), one gets a transition either in the directed 
percolation, or in the parity conserving class, depending on the choice of parameters. 

PACS number(s); 05.70.Ln, 05.70.Jk, 64.60.Ht, 02.60.Dc 



I. INTRODUCTION 

The study of systems out of thermal equilibrium has 
attracted great attention in recent years. As their equilib- 
rium counterparts, these systems may display continuous 
phase transitions characterized by a set of critical expo- 
nents Q. In particular, much interest exists in transi- 
tions from a fluctuating active state towards an absorbing 
state, i.e. a configuration where the dynamics is frozen. 

A prototype of a ID lattice model with a transition 
into an absorbing state is the contact process in 
which each lattice site can be either empty (0) or occu- 
pied by a particle {A) with the reactions A — > 2A, A ^ %. 
Depending on the relative rates of these processes, the 
stationary state is empty (when A — > dominates), or is 
occupied by a finite density of particles (if the reaction 
A — > 2 A dominates). The contact process therefore dis- 
plays a non-equilibrium phase transition which is known 
to belong to the directed percolation (DP) universality 
class. The empty lattice (00 ... 0) is the absorbing state. 

A wide range of models with transitions into absorbing 
states was found to belong to the DP universality class. 
As typical examples we quote the branching-annihilating 
random walks (BARW) with an odd number of offsprings 
1^, the Domany-Kinzel model Q and the pair contact 
process The DP class therefore appears to be ex- 
tremely robust and quite common, but it is certainly not 
the only possible one. 

Another universality class that has by now been firmly 
established is the so-called parity conserving (PC) class 
101 . A prototype model in this class is the BARW with an 
even number of offsprings Q in which particles can dif- 
fuse and undergo the reactions 2A ^ ^, A ^ [m + \)A, 
with m an even integer. In that model, the particle con- 



servation modulo two is believed to be the reason for 
the system to show non-DP critical behavior. More re- 
cently it became clear that the parity conservation, at 
least at the microscopic level, is not a necessary condi- 
tion for a PC transition to occur |7|j^. Hinrichsen |^ 
provided an example of this by introducing a one dimen- 
sional model where each lattice site can be occupied by 
at most one particle or can be in any of n inactive states 
(01, 02 ... 0n). The reactions are: 



with rate A/n (1) 

with rate /ik (2) 

with rate 1 (3) 

{k ^ I) with rate 1 (4) 



AA A%k,^kA 
A0fc,0fcyl -> 0fe0fc 

A%k,%kA~^AA 
0fc0i ^ A%i, 0feA 



We refer to the model defined by the reactions (|lj-^ as 
to the generalized contact process (GCP). The original 
contact process ||^, corresponds to the case n = 1, in 
which the reaction (||) is obviously absent. Notice that 
the reaction (0) in the case n > 2 ensures that configu- 
rations as (0i0i . . . 0i0i0j0j . . . 0j0j), with i ^ j are not 
absorbing. Such configurations do evolve in time until 
the different domains coarsen and one of the n absorb- 
ing states (0101... 0i), (0202... 02), . . . (0„0„ . . . 0,0 is 
reached. 

For the GCP with rt = 2, it was found from simulations 
1^ that the transition falls in the PC class if /ii — /i2, 
while if the symmetry between the two absorbing states 
was broken (^1 ^ ^2) & DP transition was recovered [||. 

One aim of this paper is to investigate the case n > 3, 
which has not been studied so far. Our results are ob- 
tained by a numerical investigation based on density ma- 
trix renormalization group (DMRG) techniques Q for 
the case n = 3 and n — A and an exact solution in the 
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limit — > oo. On the basis of these results we are led 
to conjecture that for n > 3 the model is always in the 
same universality class, which coincides with that of mul- 
tispecies branching and annihilating random walks ||l0| . 

There are several reasons which make the GCP inter- 
esting. Firstly, it contains the two mayor universality 
classes for transitions into absorbing states, namely DP 
(n = 1) and PC {n — 2). Secondly, for n > 1 is it inter- 
esting to study the effect of breaking the permutational 
symmetry of the model and for n > 2 the symmetry 
can be broken in different ways (see below). Finally, it 
is also interesting to investigate the performance of the 
DMRG algorithm for a system with several absorbing 
configurations and with several states per site, a situa- 
tion which is definitely more complicated than what was 
considered so far ||ll|,|2|. We recall that the application 
of the DMRG to non-equilibrium processes is rather re- 
cent and its power/limitations have not been fully inves- 
tigated yet, therefore the GCP represents another impor- 
tant testing ground for this purpose. 

This paper is organized as follows: in Sec. || we present 
the numerical results for the cases n = 1,2,3 and 4. In 
Sec. Ill we show that in the limit n oo the dynamics of 
the model can be mapped onto that of the zero tempera- 
ture n-state Potts model, which allows the determination 
of one critical exponent. Next, we present in Sec. |^ a 
conjecture for the critical behavior of the model for ji > 3 
based on the numerical and exact results. Finally, in Sec. 
^ we analyse the effect of breaking the symmetry in the 
n = 3 case, while Sec. VI concludes the paper. 



II. DMRG STUDY OF THE MODEL 

As a starting point for our analysis we use the quan- 
tum form of the master equation to describe the 
evolution of the stochastic processes in continuous time: 



d\P{t)) 
dt 



-H\P{t)) 



(5) 



where \P{t)) is a state vector whose elements are the 
probabilities of finding the system in a certain configura- 
tion, and the entries of the matrix H are the transition 
probabilities per unit of time between different configu- 
rations. As in quantum mechanics we will call H the 
Hamiltonian of the system. However, in the most gen- 
eral case, as in the GCP, the matrix H is non-hermitian, 
therefore one should distinguish between right and left 
eigenvectors, which are now not related by transposition. 
As a second consequence, the eigenvalues could be com- 
plex, but H always has at least one eigenvalue that is 
zero, and the real part of the non vanishing eigenvalues 
is strictly positive. Since we are interested in the station- 
ary behavior of the system and the relaxation towards it, 
we will determine the low lying spectrum of the Hamilto- 
nian, i.e. the part of the spectrum with the smallest real 
part. In sections II-IV we will consider the system to be 
symmetric in the ground states: iii = fi2 = ■ ■ ■ = 1^- 



First of all, the conservation of probability always en- 
sures the existence of a trivial left eigenvector with zero 
eigenvalue {0| = J2a(^\' where the sum is extended over 
all possible configurations Besides this left ground 
state, the GCP has n trivial right ground states: the 
n absorbing configurations \tp^) = |0/c0fc • ■ ■ 0fc), with 
fc = 1, 2 . . . n. To study the rest of the low lying spectrum 
we apply the DMRG in combination with finite size scal- 
ing. The DMRG is an iterative algorithm through which 
one constructs approximate eigenvalues and eigenvectors 
of H for long chains. At each iteration the lattice size 
is increased and the configurational space is truncated 
efficiently, so that one considers, instead of the exact op- 
erator iJ, effective matrices of reduced dimensions, which 
can be handled numerically. The accuracy obtained for 
the dominant eigenvalues and eigenvectors is often very 
good. Although originally invented for hermitian prob- 
lems, the DMRG also works in non-hcrmitian cases, as it 
has been shown in several examples [ pTjjT^ . 



A. The cases n = 1 and n — 2 

For n > 1 the GCP has more than one absorbing 
state and therefore more than two states per site. This, 
together with the fact that the Hamiltonian is non- 
hermitian, makes a DMRG study of the model techni- 
cally difficult. Consequently it is of great importance to 
have a way to check the results of our method. Until now 
the GCP has only been studied by means of Monte Carlo 
simulations, for n = 1,2 ||]. Therefore, this subsection 
is restricted to these two cases. We will explain how we 
applied the DMRG, give some details on the finite size 
scaling for n = 2 and compare our results with . 

A first quantity which can be calculated by the DMRG 
is the gap F between (the real part of) the eigenvalue of 
the first excited state of the Hamiltonian H and the ab- 
sorbing ground states jV'fe)- T is the inverse relaxation 
time and enables us to determine the critical region and 
the dynamical critical exponent z — v^^/v^. A direct im- 
plementation of the DMRG is however very unpractical 
since for n > 1 we have a degenerate ground state, and 
the DMRG is known to work best for gapped systems. 
We alleviated this problem by adding at the two bound- 
ary sites the following reactions 



^1 (fc 1) with rate p 



(6) 



(Recall that it is customary to use open boundary con- 
ditions in DMRG ||). With (|), only \tP^) is left as a 
ground state. The bulk critical behavior is however ex- 
pected to be unchanged. Next we performed the trans- 
formation [nsj: 



H'iA)^H + A\ri){0\ 



(7) 



with A > and where H is constructed from the reac- 
tions (0^) and (^). H' is no longer a stochastic Hamil- 
tonian, the zero eigenvalue is shifted to A, but the rest 
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of its spectrum is exactly the same as that of H. This 
can easily be checked by looking at the eigenvalues of 
the left eigenvectors which are the same for both matri- 
ces. Therefore the calculation of the gap of the original 
Hamiltonian with an n-times degenerate ground state is 
reduced to the calculation of the lowest eigenvalue of H' 
(obviously provided A is bigger than the gap of H). This 
strategy could be generalized to other systems with sev- 
eral absorbing states, provided one finds an appropriate 
boundary reaction which "eliminates" zero eigenstates, 
as done by (i) for the GCP. 
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FIG. 1. Phase diagram of the generalized contact process 
for Hi = fj,2 = ■ ■ ■ iJ-n = fJ,, as obtained by non-hermitian 
DMRG techniques, for n — 1 (filled circles), n = 2 (empty 
squares) and n = 3 (thick solid line). For n = 1 and n — 2, 
the region above the curves is the active phase, while the re- 
gion below is inactive. In the limit A ^ oo the critical line 
for n = 2 approaches a finite value of fi, while that for n = 1 
diverges towards /i — > 0. In the other limit A ^ both Unes 
merge into a common special point (see Ref. [^). For n = 3, 
the model is active in the whole parameter space except along 
the critical line = 0. 

Before turning to the finite size analysis of the calcu- 
lated gap, let us first present the phase diagram for n — 1 
and n = 2. Figure |^ shows these diagrams, obtained by 
the DMRG method and standard finite size scaling tech- 
niques. The region above the lines denotes the active 
state, where the system has a finite stationary density of 
particles, while below these lines the stationary density is 
zero and the system is in one of the absorbing states. We 
find (see below) that the critical exponents are the same 
all along the critical line and are those of the DP and 
PC universality class for n = 1 and n = 2 respectively. 
We notice that the active region increases from n = 1 to 
n = 2. The location of the critical lines agrees well with 
Monte Carlo simulations by Hinrichsen 

It is worth while at this point to present in some more 
detail the finite size scaling analysis employed for n = 2, 



in order to clarify the differences with the higher n case. 
For each choice of the rates (A, /i) we can calculate the 
gap Pi of the matrix H for a system of length L. This gap 
equals the lowest eigenvalue of H' defined in Eq. (^ . As 
a function of L, P^ has a different scaling behavior in the 
three regions of the phase diagram. In the active phase 
the gap should scale as Pi ^ exp(— aL), since asymp- 
totically in L the absorbing states are degenerate with a 
state having a finite density of particles. At the critical 
line we have P^ ~ with z = i'\\/i'± the dynamical 

exponent. For n — 2 the system is in the PC univer- 
sality class and the whole inactive phase is known to be 
algebraic: the gap decays as Pl ~ L~^. Notice that the 
latter behavior is different in DP models where the gap 
is finite, i.e. P^ ^ Fq > 0, in the inactive phase. 

The finite size scaling behavior of P^ can best be mon- 
itored by plotting the discrete logarithmic derivative of 
the gap: Yl = ln(Pi+i/Pi_i)/ ln[(L + 1)/(L - I)]. In 
the active phase the scaling form P^ ~ exp (—aL) implies 
Yl —aL (with a a positive constant), while Yl ^ —z 
at the critical line. 
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FIG. 2. Plot of Yl for the GCP with n = 2 as function of 
fi and for A = 0.5. Inset: plot of zl as function of the inverse 
system length 1/L. The horizontal lines indicate the values 2 
and the known value of the exponent z for the PC universality 
class (2 = 1.74). 

Figure || shows a plot of Yl for n = 2,A = 0.5as 
a function of the parameter fi for L = 5, 7 ... 23 (in the 
present case the calculation of the gap was extended up to 
L = 24). In the right part of the figure one distinguishes 
the scaling for the active phase, while in the left part 
(large /i) one sees that Yl flattens out and approaches 
the value of —2, as expected for the algebraic inactive 
phase of PC models. The maxima of Yl, marked with 
a filled circle in Fig. |^, identify the boundary between 
the active and inactive phase and can be used as critical 
point estimates. From an L — > 00 extrapolation we find 
/ic = 0.69(1) for our choice of A = 0.5, which determines 
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point B of the critical line in Fig. |l|. In the inset of Fig. ^ 
we plot zl = — max^ A = 0.5), as function of 1/L. 

Extrapolating in the limit L — > oo we find zl 1.14,1 ^ 
in accurate agreement with the corresponding exponent 
known for the PC class z = vw/vi. = 1.74(1) |16[ . 

In this way, the DMRG enables us to determine the 
critical region and the exponent z. As we will show in 
the following subsection, using different boundary condi- 
tions we can also find estimates for other exponents. In 
Table | we show our estimates of the critical exponents 
z and P/vx, in the points A and B of Fig. ^ We re- 
cover the exponents of DP and PC as expected for n = 1 
and n = 2 respectively, indicating that the DMRG is 
performing well. 

TABLE I. Estimates for the critical exponents z = 
and l3/v±_ calculated for n = 1 and 2 absorbing states. We re- 
call that for DP z = 1.5806, l3/ux_ = 0.25 and for PC z = 1.75, 
l3/vi_ = 0.50. 



n 


A 


Mc 


z = i^ii 




1 


0.5 


0.20(1) 


1.575(5) 


0.255(5) 


2 


0.5 


0.69(1) 


1.747(5) 


0.49(1) 



B. The case n = 3 

Since the results of the DMRG are consistent with the 
Monte Carlo results for n = 1 and n = 2, we can confi- 
dently use it to study also the case n = 3. Again we start 
by an analysis of the gap: in Fig. ^ we plot the quan- 
tity Yl along the hue A = for L = 7, 9 . . . 19. As for 
n = 2 we find a clear evidence of an active phase for small 
fi where —aL. However, when we want to deter- 

mine the critical point by localising the maxima, we find 
that upon increasing L the maxima shift towards larger 
values of /i. In the limit L — > cxd we find //max ~* oo, 
indicating that the system is always active and that the 
critical point is shifted to l//i = 0. This numerical evi- 
dence together with the arguments presented in the next 
section lead us to expect that the model is critical only 
if l//i = 0, as indicated in Fig. |l|. The estimate of the 
exponent z at X = n —>■ oo yields z = 2.00(3). In the 
next section we will give an analytical treatment of the 
case l//i = confirming this value of z. 

To get access to more critical exponents, we introduce 
different boundary conditions. We replace (||) with: 



A with rate p' V k 



(8) 



through which particles are continuously injected at the 
boundary sites. Now the states \tp'^) are no longer absorb- 
ing configurations, and the spectrum of the Hamiltonian 
is non-degenerate. There is a unique right eigenvector 
\(t>o)i which for finite system lengths always has a finite 
density of particles. In the active phase the gap F^ is 
no longer asymptotically degenerate, as in the previous 



case, but it remains finite, indicating that the system 
relaxes exponentially fast towards |0o)- Summarizing, 
if we consider the model defined by the reactions (|l[]^) 
and the boundary term (^) we have a gap behaving as 
limL^oo F'^ = Fq > in the active phase, while as before 
F^ ~ at the critical point. This means we can use 
F^ itself to discriminate between the active and critical 
domains. 
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FIG. 3. Plot of Yl for the GCP with n = 3 along the line 
X = fi. 




FIG. 4. Plot of the gap r'£(/t), calculated along the line 



for 



\ = ji (with p' = /i in reaction (^)) as function of ^ ^ 
systems of lengths L = 6, 8 ... 18. Inset: scaling collapse of 
/i''liri,(Ai), plotted as function of /iL"^/""^ for L = 10, 12, 14, 
16 and 18 and where we used v« = 2 and iy± = 1. 



Figure ^ shows a plot of the gap F^ along the line 
X — fi. When we extrapolate L — > oo, the gap remains 
finite for every value of /x, i.e. we find again that the 
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system is active throughout the whole parameter space. 

These data also give the possibility to estimate the 
correlation length exponent along the time direction: v\\ . 
Finite size scaling around the critical point 1//^ = gives 
us the following scaling relation: 



Pc«{^;L) = - 



(9) 



For fixed and L oo, remains finite, so the scal- 
ing function should behave like lim2,^g+ f{x) = fo- As a 
consequence in the thermodynamic limit L oo the gap 
should vanish at the critical point as F^ ~ . The 
linear behavior of the gap in Fig. ^, when plotted as a 
fimction of already indicates that = 2 and this re- 
sults is also consistent with the scaling collapse reported 
in the inset, where we took = 2 and z = 2. 

Finally using the same boundary term (H) it is also pos- 
sible to compute the critical exponent P which describes 
the behavior of the particle density near the critical point. 
With the boundary condition (||) the stationary state of 
any finite system has a non-zero density of particles. We 
calculated in particular the average density of particles 
in the central site of the chain which we will denote as 
PL(A,/i). This quantity, calculated along the line A = /i, 
is shown in Fig. |^ for chains of various lengths (up to 
L = 24) and plotted as function of 1/^. For large fj, and 
A (/i > 2) the DMRG does not perform so well and we 
had to restrict the calculation to L = 14. 
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FIG. 5. Plot of the stationary particle density p for the 
model with n = 3 absorbing states along the line fi = X for 
various system sizes (L = 4, 6, . . . 24). For large fj, the calcu- 
lation has been limited to L = 14. Inset: plot of the effective 
exponent /3 as function of l//i for L — 6 (lower curve) up to 
L — 14 (upper curve). 

As seen in the figure in the limit L — > oo the sta- 
tionary particle density vanishes only for fi oo, again 
indicating that the system is always active for any finite 
/i. For large /i the order parameter is expected to decay 
as p ^ therefore the quantity: 



d\npL{\ ^ p) 
din p 



(10) 



converges to /3 for L — > oo and /i — > oo. The inset of 
Fig. H shows a plot of PcS as function of 1 //i for chains 
of various system lengths. The extrapolated result for 
L = 12 and L = 14 is /3 = 1.00(1). 

In summary, we found that for n — i, the system is 
always active, except for l//i = where it is critical and 
we have exponents consistent with z = 2, v\\ = 2, [3 — 1. 



C. The case n = 4 

To conclude our numerical calculations for the case of 
symmetrical ground states, we studied the system with 
four inactive states per site. Here we restricted ourselves 
to a calculation of the density pi(A,/i) with boundary 
condition (0) . Our results are shown in figure 0. 




FIG. 6. Plot of the stationary particle density p for the 
model with n = 4 absorbing states along the line p. = X for 
various system sizes (L = 4, 6, . . . 14). For large p the cal- 
culation has been extended up to L = 10. Inset: plot of the 
effective exponent (3 as function of 1 /p for L = 6 (lower curve) 
up to L = 10 (upper curve). 

From these data we find, as was the case for n = 3 
that the system is active for any finite value of p, and 
that the critical exponent /? = 0.98(4). 



III. FAST RATE EXPANSION FOR OR A ^ oo 



In this section we show how the dynamics of the model 
can be simplified in the limit where one of the rates of the 
decay processes (^ or (||) goes to infinity. For p ^ oo 
the resulting effective dynamics coincides with that of a 
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zero temperature n-state Potts model. This leads to a 
determination of the exponent z. 

Intuitively, the argument goes as follows. For /i — s- cx), 
all particles in the system will disappear very quickly, and 
in any finite system one will soon be in a configuration 
containing only the n empty sites 0i, 02 and 0„. Particles 
will then be created at the boundaries separating inactive 
domains by (^, but again they will disappear immedi- 
ately by reaction (|^). Hence, when looked at the time 
scale of the slow processes (^) and (||) the dynamics can 
be limited to the set of configurations that contain only 
empty sites. (In this limit the creation of two particles 
on nearest neighbour sites occurs with even smaller prob- 
ability and the effect of can therefore be neglected). 
As an example consider: 



(4) 



''kVkVkVWl ■ ■ 



(2) 



UkUkUlVlVl . ■ 



(11) 



The whole process consists of moving the domain wall 
one unit to the left with rate (on the time scale of the 
slow process) equal to 1/2. Going in this way through 
all possibilities one can derive an effective dynamics on 
the slow time scale, which turns out to involve diffusion, 
annihilation and coagulation of domain walls. 

The above heuristic reasoning can be made mathemat- 
ically rigorous using a fast rate expansion introduced in 
[llil . We therefore write H as 



H = Ho + nHi 



(12) 



where Hq contains the reactions (|l|) , (|^) and (^) while Hi 
contains only the reactions (|^). In the limit /i — > oo, it 
is appropriate to apply the Schwinger-Dyson formula to 
the operator e~^* which appears in the formal solution 
of the master equation (||). One has 



1- / dTiHoin) 



dTl f ' dT2Ho{Tl)Ho{T2) 

Jo 



where 



Ho{t) = e^^'^'^Hoer""'^ 



(13) 



(14) 



Using this expansion it can be shown that for /i — > oo 
the time evolution operator e"''^' becomes 



lim e-(^°+^^i)* = e-ffoty* 

fl — >-oo 



(15) 



where 

Ho = T*HoT* (16) 
and T* is the projector on the ground states of Hi, i.e. 

(17) 



T" = lim e-"^* 



Equation ( p^ shows that for /i ^ oo the generator of 
the effective dynamics is Hq which is nothing but Hq 
projected onto the ground states of Hi. In the GCP, 
the ground states of Hi for a chain of length L are the 
configurations without particles and Hq is then the 
effective Hamiltonian of the slow processes projected in 
this reduced space. This is the mathematical description 
of the physical arguments given in the beginning of this 
section. If one works out the matrix elements of Hq (see 
the appendix) one obtains that the following processes 
can occur with the indicated rates 



''kHIlVk 
%kUl 



'ikVlVm 



hMk rate 2 
h%k%i rate 1/2 
rate 1/2 
h^k^m rate 1 
h^r,Sm rate 1 



(18) 
(19) 
(20) 



It is now useful to interpret the n empty states as the 
possible spin values of an n-state Potts model. In this 
language, the processes ( p"8[ - |20| ) can be summarized as 
follows: the central spin assumes the value of one of its 
neighbour spins with equal probability. The dynamics 
of our model in the limit /i — > oo is therefore consistent 
with the requirements of detailed balance for an n-state 
Potts model at zero temperature. It is generally expected 
that if such a dynamics includes a domain wall diffusion, 
that then it is critical with a dynamic exponent z — 2 
|]l7| , independently of n. In fact, the exponent z = 2 
can be derived exactly for the case that the rate of the 
process (|l8) equals 1, and those of the processes in (|20| ) 
equal 1/2 |18[. Hence, it is quite possible that also for 
our model, z = 2 exactly. For n = 2 we thus recover the 
known dynamical exponent z = 2 in the inactive phase of 
a model with a PC transition. For n > 3 our numerical 
data strongly suggest that for /i ^ oo the GCP is crit- 
ical, and hence the exponent z = 2 must correspond to 
the dynamical exponent at criticality for these models. 
This estimate is indeed consistent with the value that 
was determined numerically for n = 3 in the previous 
section. 

The result that numerically the exponent (3 is the same 
for n = 3 and n = 4, combined with the fact that z = 2 if 
n > 3 leads us to conjecture that for all n > 3 the critical 
exponents are the same, i.e. /3 = 1, z = 2 and = 2. 
In the next section we will give further arguments that 
support this conjecture. 

It is also possible to study the effective dynamics of 
the GCP for A ^ oo. In this limit, the dynamics of 
the model when considered on the time scale of the slow 
processes (^||), will be restricted to the space of config- 
urations without particle pairs. Each particle present in 
the system then separates two domains of empty sites. 
Therefore particles can be labelled by two indices and in 
this way they can be divided into classes. We will de- 
note by Aij a particle A that separates an 0i-domain on 
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its left side from an 0j-domain on its right side. In the reactions for n = 2 are: 
limit A ^ oo, it is useful to look at the dynamics of these 
particles. When n = 1 there is only one type of particle, 
and from a determination of the effective Hamiltonian _ffo 
(see the appendix) we find that this particle can diffuse, 
and undergo the reactions 2 A A and A ^ %. Since 
there are no processes that create particles, we arrive at 
the conclusion that independently of /i, the GCP with 
n = 1 must always be in the inactive state when \ ^ oo. 
This is in agreement with our numerical results (see Fig. 

For n > 1, the situation is less clear. Now there are 
processes that both destroy and create particles present 
in the effective dynamics (see the appendix). It is there- 
fore in principle possible to have both an active and an 
inactive phase, depending on the value of /i. At this mo- 
ment, we can draw no firm conclusions for the form of 
the phase diagram when A — > oo. On the basis of our nu- 
merical work and on the basis of our results for /i — > cxd, 
we believe that the model is probably always active along 
that line, at least when n > 2. 



(21) 



(for n — 2, one has obviously i = 1, j = 2 or vicev- 
ersa). An example of such reactions is shown in Fig. |^ 
(a) and (b). Notice that the reactions (a) and (b) given in 
( pi| ) are those for a branching-annihilating random walk 
(BARW) with two offsprings, which suggests indeed, as 
found numerically that the universality class is PC. 

These arguments can be extended to the case rt > 2, 
where there is still the possibility of having reactions of 
type (a) and (b), but also reactions involving three dif- 
ferent domains {i ^ j, i ^ k and j ^ k).: 



Xr 



XikXkj 



XihX 



(d) 



Xj, 



(22) 





(b) 




(d) 



space 

> 



When n > 2 there are actually n{n — l)/2 domain walls, 
and the model described by reactions (|2^) and ( p^ ) is 
now a BARW with more than one type of particles. To 
our knowledge this type of model has not been studied 
yet, but we expect it to be in the same universality class 
as the GCP with n > 2, i.e. always in an active state 
except when the rates for the processes (a) and (c) are 
zero, and with exponents z — 2, f3 = I and = 2. 

There is a BARW with more than one type of particles 
that has attracted some attention recently. The model 
was introduced by Cardy and Tauber [nol, who consid- 
ered a system with N different particles A°', where a = 1, 
2, . . . N which diffuse and undergo the reactions: 



time 



2A" 

A" ^ A°'A°'A" 
A" ^ A°'A'^A'^ 



with rate 1 
with rate a 
with rate ct7(N 



(23) 
(24) 
(25) 



FIG. 7. Possible reactions in the coarse-grained represen- 
tation of the GCP for n — 2 (a-b) and n > 2 (a-b-c-d). 



IV. RELATION WITH BRANCHING AND 
ANNIHILATING WALKS 

In his paper, Hinrichsen |^ gave an heuristic argument 
which relates the GCP for n = 2 to the BARW with two 
offsprings, thus explaining the PC universality found nu- 
merically. This argument works as follows: indicating 
with Xij the domain wall between two configurations 0; 
and 0j he considered an effective dynamics for the vari- 
ables Xij . As seen in the previous section such represen- 
tation of the original model becomes exact either in the 
limit X ^ oo and then Xij coincides with the particle 
Aij , or in the limit /i — > oo where Xij coincides with the 
bond variable 0^0^. For finite values of these parame- 
ters one can still apply this reasoning at a coarse-grained 
level. In this case Xij is not a sharp domain wall, but an 
object with a fluctuating thickness. Hinrichsen argued 
that in this coarse-grained representation the most likely 



where in the last reaction it is understood that a (3. 
The bosonic version of this model has precisely the same 
exponents we determined for the GCP with n > 2 [p^ . 
Notice that the coarse-grained representation of the GCP 
as defined by the reactions in Fig. ^ and the model de- 
fined by the reactions (|2^j2^ , p5| ) do not actually coincide. 
While there is an obvious correspondence between ([2^), 
(H) with (b), (a) of Fig. 0, the reaction (|5|) does not 
have any obvious counterpart. It is not a priori clear 
therefore that the two models are in the same universality 
class. The coincidence of the critical exponents therefore 
suggests that one could replace (|5|) with other reactions, 
as for instance A" A'^A^, without changing the uni- 
versality class. To our knowlegde, BARW models with 
this kind of reactions have not been studied yet. They 
form an interesting subject for further investigation. 

The model of Cardy and Tauber has raised some in- 
terest recently since it has been found that fermionic and 
bosonic versions of the model are in different universality 
classes [ p^ . In the fermionic version only one particle per 
lattice site is allowed, which implies that particles of dif- 
ferent species block each other. In the fermionic version 
of the model it makes for instance a difference wether the 
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two offsprings produced by the reaction (g5|) are placed 
to the same side or at opposite sides of the parent par- 
ticle II^JI^. If, for instance, they are placed at opposite 
sides the offsprings cannot annihilate through the reac- 
tion ( p3| ) because of the presence of the parent particle 
that blocks them. In the bosonic model the two offspring 
can instead always recombine. 

It is important to stress here that in the multispecies 
BARW model we constructed ((^-(p2|)) from a coarse- 
grained representation of the GCP there are no blocking 
effects. By the very construction of the model two do- 
main walls approaching each other can always annihilate. 
Therefore even if the model is clearly of fermionic nature 
its universality class, as found numerically, is that of the 
bosonic multispecies BARW. 



V. THE EFFECT OF BREAKING THE 
SYMMETRY 

As a final point we consider the effect of breaking 
the permutation symmetry of the inactive states of the 
model. For n = 2 Hinrichsen Q explicitely broke the Z2 
symmetry by choosing /ii 7^ in the reaction (H). As 
a result he found the system to switch from PC to DP 
behavior, which was understood as PC being related to 
the presence of an exact Z2 symmetry. 



1.5 




0.5 




PC 



0.05 0.1 0.15 0.2 0.25 
1/L 



0.05 0.1 0.15 0.2 0.25 
1/L 



FIG. 8. Plot of 5 for the GCP with n = 3 and (a) 
/Lti/2 = /i2 = M3) (b) = /i2 = Ms- Both cases were cal- 
culated with A = 0.5 and p' = 1.5. Symbols refer to curves 
calculated in the inactive phase (stars), at the critical point 
(filled circles) and in the active phase (empty squares). At 
criticality converges towards the values oi [3/v±^ expected 
for DP (a) and PC (b), and indicated by horizontal dashed 
lines in the figure. Notice the two distinct behaviors of 5 in 
the inactive phase. 

For n = 3, we now perform a similar symmetry- 
breaking by considering the following two cases: (a) 



Ml/2 = /Lt2 = M3 and (b) 2^1 = M2 = M3- In both cases 
the system has a Z2 symmetry due to the equivalence 
of the states 02 and 03. The difference is that starting 
from a random configuration in the first case the system 
is more likely to reach the absorbing states (02, 02, 02, ■ • ■) 
and (03,03,03,...) compared with (0i, 0i, 0i, . . .), while 
in the second case the reverse is true. We calculated 
the particle density pl(A, ^i, M2j Ms) before, using the 



boundary-term (f 



Depending on the phase in which 
the model is, the density will behave as = po + Ce~°^ 
(in the active region), ^ L^^/"^ (at criticality), 
or PL ^ e~°-^ (in the inactive phase). If we define 
S{L) = - In [pl+i/pl] I In {{L + 1)/ ln(L - 1)], we expect 
limL^oo <5(L) to be zero in the active phase, to be /3/i^_l 
at critical points, and -l-oo in the inactive phase. In Fig. 
^ we plotted (5(L) as a function of 1/L for cases (a) and 
(b) with the choice A = 0.5 and injection rate p' — 1.5. 
In both cases for small p one finds the typical scaling 
behavior of the active phase with b{V) — > just as for 
the case pi = P2 = fJ-a- For large pk however the situa- 
tion differs from the symmetric model. For (a), we find 
6(L) +00 for large p^ i.e. one has a standard inactive 
phase with a particle density exponentially small in L. 
In case (b), S{L) becomes equal to 1, implying that the 
inactive phase is itself critical with /3/i^j_ = 1. In between 
the active and the inactive phase we have a critical point 
where S{L) is going to a distinct finite value. The critical 
point estimates of 6{L) are marked by filled circles in Fig. 
^. For case (a) the critical point is at pi « 0.64 , while 
for (b) it is at pi ~ 0.42. The value of (i/v± agrees with 
that of DP and PC respectively, as can be seen in Fig. ^ 
where the critical values of these universality classes are 
indicated with a dotted line. 

This indicates that on breaking the symmetry of the 
inactive states the remaining symmetry of the dominant 
rates pi determines the critical behavior. In case (b) 
Ml < ^2 = Ms fli6 dominant rates still have a Z2 symme- 
try leading to the PC universality class, while in case (a) 
Ml > /^2 = M3) -DP critical behavior is recovered. 



VI. CONCLUSIONS 

In this paper, we studied a generalised contact pro- 
cess first introduced by Hinrichsen. The major part of 
our results were obtained by applying the DMRG to the 
model. With this technique we verified that for n = 2, 
the critical line of the model is in the PC universality 
class, consistent with earlier results coming from simula- 
tions. 

A first set of new results was obtained for the case 
n = 3, for which we found that the model is always active, 
except when p — > 00, which corresponds to the critical 
line of the model. From our numerical work we deter- 
mined the critical exponents to be equal to z = 2, v^^ — 2 
and (3—1. Using well estabhshed scaling laws pl| , 
other exponents can be determined from these three. For 
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n = 4, we found evidence that the phase diagram is the 
same, and that the critical exponent /? also equals 1. Sec- 
ondly, using a fast rate expansion that becomes exact for 
/i — > oo, we were able to argue that in that limit z = 2. 
It can be hoped that by examining the model for /i^^ 
small using perturbation techniques, it may be possible 
to determine also the exponents (3 and i^y exactly. On the 
basis of these numerical and exact results we conjectured 
that the universality class of the model is the same for 
all n > 3. 

The exponent values that we found for n > 3 coincide 
with those of the BARW model with more than one type 
of particles introduced by Cardy and Tauber [0. We 
were able to give an heuristic argument that explains why 
the two models could be in the same universality class. It 
is interesting to remark that despite many attempts the 
number of universality classes found for phase transitions 
out of an adsorbing phase, remains very limited. It could 
have been hoped a priori that for the generalised GCP 
studied here, new universality classes could appear for 
n > 2. In a sense our results show that this is true, but 
only in the least exciting way possible: the universality 
class does not depend on n, and moreover the exponents 
take on rather trivial values. One could hope that by 
lowering the permutation symmetry to a .Z„-symmetry, 
other universality classes could appear for n > 4. This 
could be done, e.g. by having the rates of the process 
depend on |A: — Z| mod(n). But since this would only 
make a difference for n > 4, it will probably be difficult to 
investigate such a model with the numerical techniques 
currently available. 

We also verified that if one breaks the permutation 
symmetry of the model with n = 3, one recovers a DP 
or PC universality, suggesting that it is the symmetry of 
the largest rates that determines the universality class. 

Finally, we remark that the consistency of the DMRG 
results with those coming from simulation for n — 2, or 
with the exactly determined value of z for n > 3, shows 
convincingly that the DMRG can be trusted as a power- 
ful method in the study of (criticality in) non equilibrium 
systems. 



ACKNOWLEDGMENTS 



JH and CV thank the lUAP, Belgium for financial 
support. EC is grateful to INFM for financial support 
through PAIS 1999. 



APPENDIX 

In this appendix we will calculate explicitly 
linv^oo e^*-^"^''^^-** where r is one of the rates of the 
GCP. In section III it was shown that this reduces to 
the calculation of Hq = T*HoT*, the Hamiltion that de- 



scribes the effective dynamics. Since T* = limt^oo e" 
this is a projection of Hq on the ground states of Hi. 

Let us denote the right ground states of Hi as and 
the left ones as {pi\. 



Hi It/.,) = 

{p^\Hl=0 



(26) 
(27) 



The physical meaning of \tpi) is evident: they are the 
stationary states of Hi, any state {ip) will under the dy- 
namics of Hi relax into one of these . The left ground 
states {pi I can be interpreted as linear functionals giving 
the corresponding probabilities: any state \^) will under 
the dynamics of Hi relax into the ground state \4>i) with 
probability {pi \ip) (where we assumed that {pi\ are nor- 
malized: (pilipj) = Sij). Using this notation we can write 
the projection operator T* as 



(28) 



Since clearly [^,. {pi\] \tp) = 1 for any (normalized) state 
IV'), this projection conserves probability, meaning that 
when Hq is a stochastic operator, so is Hq. Furthermore 
we can write the transition rates of the effective dynamics 
between state \4>i) and as 



rate(|V'«)^|V'j)) = (Pjl(-i?o)|^. 



(29) 



These matrix elements determine the effective dynamics, 
and we will now calculate them explicitly. 

We first study the limit where the rate /i of the pro- 
cesses j40fc,0fe^ 9>k^k goes to infinity. The Hamilto- 
nian is of the form: 



H = Ho+ fiHi 



(30) 



where Hq is the generator of the processes (|i|), (||), (^) 
and Hi is the generator of process (^) with the factor fi 
brought out. In this case the ground states of Hi 
are all configurations containing no particles A. Firstly, 
we note that the processes (|l|) and (||) can only act on 
configurations containing particles, so they cannot con- 
tribute to the rates ( p^ ) of the effective dynamics and we 
redefine Hq without them. 

When we now project Hq which contains only 2-site 
interactions, onto the the resulting operator will 

contain 3-site interactions. It is therefore convenient to 
first rewrite Hq as a 3-site operator. Since we only have 
reaction (^) left in Hq, this becomes {k ^ I ^ m) 



OkVWk 



> ^kA^k rate 2 
%kA^rn rate 2 
■ 0fe^0/ rate 1 
%kAA rate 1 



(31) 



(together with some reactions that are obtained by re- 
flection). We finally notice that the last process of ( |3l| ) 
is again not relevant for the projection on |'0i)j and de- 
termine the effective dynamics in the following diagram 
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reaction 
with rate.. 



projection 
with probabihty.. 



netto rate 







1 ^ 


0fc( 




rate 2.1 


= 2 












rate 2.^ 


= 1 (32) 












rate 2.-i 


= 1 






1/2 


<l>k 


0fc0i 


rate 1.^ 


1 

2 






1/2 




0z0i 


rate 1.^ 


_ 1 
2 



which are the processes and their corresponding rates 
already given in (p^po|). Note that this calculation is 
exactly the same for /j, = A — > oo. 

Next, we consider now the limit A — > oo. In this case 
we have 



H = Hq + A-ffi 



(33) 



where Hq is the generator of the processes (^^, and Hi 
is the generator of process (|^) with the factor A brought 
out. The ground states of Hi are now all configurations 
containing no particle pairs. In contrast with the pre- 
vious case, all processes of Hq are now relevant for the 
projection on the ground states of Hi. 

We will start with the case n — 1, where there is only 
one inactive state 0, and process can a priori not take 
place. For process (^) we again use the 3-site represen- 
tation, while process is so simple that we keep the 
2-site representation. We then get the following effective 
dynamics 



reaction 
with rate.. 



projection 
with probability.. 



A0 



rate fj. 



AAd) 



AU AAA 



1 

1/2 
1/2 
1/2 
1/4 
1/8 
1/8 



0A0 
A%A 

%A% 

Am 

UA 



netto rate 



rate /z.l = 
rate 1-5 = ^ 
rate ^■\ — \ 
rate 2.\ = 1 



(34) 



rate 2.i 
rate 2.i 
rate 2.i 



For A — > cx) and n = \ we find the effective dynamics to 
contain only diffusion and destruction of particles. Be- 
cause of the first reaction appearing in (|34|), the decay 
of particles is exponentially fast, meaning that for any 
finite value of /i this system is non-critical. 

For the case n > 1 also process (^ has to be taken 
into account. As a consequence, two different neighbour- 
ing inactive domains remain active. For example 

MSi ^ hAd^i (35) 

remains a process of the effective dynamics. One can eas- 
ily construct the complete effective dynamics, but this 
does not lead to much further insight in the phase dia- 
gram. One can only conclude that because of the pres- 
ence of the process (|35|), it is in principle possible that 
both active and inactive phases are present. 
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